# -*- Mode: shell-script -*- 

#############################################################################
##
#A  ctfilter.g                  GAP library                     Thomas Breuer
##
##
#Y  Copyright (C) 2018-2021, Carnegie Mellon University
#Y  All rights reserved.  See LICENSE for details.
#Y  
#Y  This work is based on GAP version 3, with some files from version 4.  GAP is
#Y  Copyright (C) (1987--2021) by the GAP Group (www.gap-system.org).
##
##  This file contains those functions that are used as filters for
##  characters.
##
##


#############################################################################
##
#F  StepModGauss( <matrix>, <moduls>, <nonzerocol>, <col> )
##
##  Gauss elimination for column <col> of <matrix>, where the entries of column
##  'i' are taken modulo '<moduls>[i]' and only those columns 'i' with
##  '<nonzerocol>[i] = true' (may) have nonzero entries;
##
##  after that the only row containing a nonzero entry in column <col> will be
##  the first row of <matrix>, and again Gauss elimination is done with that
##  row and the row $\delta_{<'col'>}$;
##  if there is a row with nonzero entry in column <col> this row is returned,
##  otherwise 'false' is returned.
##
StepModGauss := function( matrix, moduls, nonzerocol, col )
    local i, k, x, y, z, a, b, c, d, val, stepmodgauss;

    if matrix = [] then return false; fi;
    matrix[1][col]:= matrix[1][col] mod moduls[col];
    for i in [ 2 .. Length( matrix ) ] do
      matrix[i][col]:= matrix[i][col] mod moduls[col];
      if matrix[i][col] <> 0 then     # transform so that matrix[i][col] = 0;
        z:= Gcdex( matrix[1][ col ], matrix[i][col] );
        a:= z.coeff1; b:= z.coeff2; c:= z.coeff3; d:= z.coeff4;
        for k in [ 1 .. Length( nonzerocol ) ] do
          if nonzerocol[k] then
            val:= matrix[1][k];
            matrix[1][k]:= ( a * val + b * matrix[i][k] ) mod moduls[k];
            matrix[i][k]:= ( c * val + d * matrix[i][k] ) mod moduls[k];
          fi;
        od;
      fi;
    od;
    if matrix[1][col] = 0 then return false; fi;  # col has only zero entries
    z:= Gcdex( matrix[1][col], moduls[col] );
    a:= z.coeff1; b:= z.coeff2; c:= z.coeff3;
    stepmodgauss:= [];
    for i in [ 1 .. Length( nonzerocol ) ] do
      if nonzerocol[i] then
        stepmodgauss[i]:= ( a * matrix[1][i] ) mod moduls[i];
        matrix[1][i]:= ( c * matrix[1][i] ) mod moduls[i];
      else
        stepmodgauss[i]:= 0;
      fi;
    od;
    stepmodgauss[col]:= z.gcd;
    matrix[1][col]:= 0;
    return stepmodgauss;
    end;


#############################################################################
##
#F  ModGauss( <matrix>, <moduls> )
##
##  <matrix> is transformed to an upper triangular matrix generating the same
##  lattice modulo that generated by
##  $\{<moduls>[i] \cdot \delta_i; 1 \leq i \leq \|<moduls>\|\}$.
##
##  <matrix> will be destructed, the triangular matrix is returned.
##
ModGauss := function( matrix, moduls )
    local i, modgauss, nonzerocol, row;
    modgauss:= [];
    nonzerocol:= [];
    for i in [ 1 .. Length( moduls ) ] do nonzerocol[i]:= true; od;
    for i in [ 1 .. Length( matrix[1] ) ] do
      row:= StepModGauss( matrix, moduls, nonzerocol, i );
      if row <> false then Add( modgauss, row ); fi;
      nonzerocol[i]:= false;
    od;
    return modgauss;
    end;


#############################################################################
##
#F  ContainedDecomposables(<constituents>,<moduls>,<paracharacter>,<func>)
##
##  <constituents> must be rational vectors, <paracharacter> a
##  parametrized rational vector. Using 'StepModGauss' all elements $\chi$ of
##  <paracharacter> are calculated which modulo <moduls> lie in the lattice
##  spanned by <constituents> and satisfy $<func>( \chi )$.
##
ContainedDecomposables := function(constituents,moduls,paracharacter,func)
    local i, x, matrix, fusion, newmoduls, candidate, classes,
          nonzerocol, possibilities, images, uniques,
          nccl, min_anzahl, min_class, erase_uniques, impossible, 
          evaluate, remain, ncha, pos, fusionperm, newimages, oldrows,
          newmatrix, step, erster, descendclass, j, row, oldimages;
    
    # step 1: check and improve input
    
    if IsList( paracharacter[1] ) then     # necessary if no class is unique
      min_anzahl:= Length( paracharacter[1] );
      min_class:= 1;
    fi;
    matrix:= CollapsedMat( constituents, [ ] );
    fusion:= matrix.fusion;
    matrix:= matrix.mat;
    newmoduls:= [];
    for i in [ 1 .. Length( fusion ) ] do
      if IsBound( newmoduls[ fusion[i] ] ) then
        newmoduls[ fusion[i] ]:= Maximum( newmoduls[ fusion[i] ], moduls[i] );
      else
        newmoduls[ fusion[i] ]:= moduls[i];
      fi;
    od;
    moduls:= newmoduls;
    candidate:= [];
    nonzerocol:= [];
    for i in [ 1 .. Length( moduls ) ] do
      candidate[i]:= 0;
      nonzerocol[i]:= true;
    od;
    possibilities:= [];  # global list of all vectors $\chi$ in the lattice
                         # which satisfy $'func( \chi ) = true'$
    images:= [];
    uniques:= [];
    for i in [ 1 .. Length( fusion ) ] do
      if IsInt( paracharacter[i] ) then
        if ( IsBound( images[ fusion[i] ] ) ) then
          if IsInt( images[ fusion[i] ] ) and
             paracharacter[i] <> images[ fusion[i] ] then
            return [];
          elif IsList( images[ fusion[i] ] ) then
            if not paracharacter[i] in images[ fusion[i] ] then
              return [];
            else
              images[ fusion[i] ]:= paracharacter[i];
              AddSet( uniques, fusion[i] );
            fi;
          fi;
        else
          images[ fusion[i] ]:= paracharacter[i];
          AddSet( uniques, fusion[i] );
        fi;
      else            # IsList( paracharacter[i] )
        if not IsBound( images[ fusion[i] ] ) then
          images[ fusion[i] ]:= paracharacter[i];
        elif IsInt( images[ fusion[i] ] ) then
          if not images[fusion[i]] in paracharacter[i] then return []; fi;
        else          # IsList
          images[ fusion[i] ]:=
                      Intersection( paracharacter[i], images[ fusion[i] ] );
          if images[ fusion[i] ] = [] then
            return [];
          elif Length( images[fusion[i]] ) = 1 then
            images[ fusion[i] ]:= images[ fusion[i] ][1];
            AddSet( uniques, fusion[i] );
          fi;
        fi;
      fi;
    od;
    nccl:= Length( moduls );
    
    # step 2: first elimination before backtrack:
    
    erase_uniques:= function( uniques, nonzerocol, candidate, images )
    # eliminate all unique columns, adapt nonzerocol;
    # then look if other columns become unique or if a contradiction occurs;
    # also look at which column the least number of values is left
    local i, j, abgespalten, col, row, quot, val, ggt, a, b, k, u,
    firstallowed, step, gencharacter, newvalues;
    abgespalten:= [];
    while uniques <> [] do
      for col in uniques do
        candidate[col]:= ( candidate[col] + images[col] ) mod moduls[col];
        row:= StepModGauss( matrix, moduls, nonzerocol, col );
        if row <> false then
          abgespalten[ Length( abgespalten ) + 1 ]:= row;
          if candidate[ col ] mod row[ col ] <> 0 then
            impossible:= true;
            return abgespalten;
          fi;
          quot:= candidate[col] / row[col];
          for j in [ 1 .. nccl ] do
            if nonzerocol[j] then
              candidate[j]:= ( candidate[j] - quot * row[j] ) mod moduls[j];
            fi;
          od;
        elif candidate[ col ] <> 0 then
          impossible:= true;
          return abgespalten;
        fi;
        nonzerocol[ col ]:= false;
      od;
      min_anzahl:= "infinity";
      uniques:= [];
      for i in [ 1 .. nccl ] do
        if nonzerocol[i] then
          val:= moduls[i];
          for j in [ 1 .. Length( matrix ) ] do # zerocol iff val = moduls[i]
            val:= GcdInt( val, matrix[j][i] );
          od;
    
    # update lists of image
    
          newvalues:= [];
          for j in images[i] do
            if ( candidate[i] + j ) mod val = 0 then
              AddSet( newvalues, j );
            fi;
          od;
          if newvalues = [] then                   # contradiction
            impossible:= true;
            return abgespalten;
          elif Length( newvalues ) = 1 then        # unique
            images[i]:= newvalues[1];
            AddSet( uniques, i );
          else
            images[i]:= newvalues;
            if Length( newvalues ) < min_anzahl then
              min_anzahl:= Length( newvalues );
              min_class:= i;
            fi;
          fi;
        fi;
      od;
    od;
    if min_anzahl = "infinity" then
      gencharacter:= Indirected( images, fusion );
      if func( gencharacter ) then Add( possibilities, gencharacter ); fi;
      impossible:= true;
    else
      impossible:= false;
    fi;
    return abgespalten;
    # impossible = true: calling function will return from backtrack
    # impossible = false: then min_class < "infinity", and images[min_class]
    #                     contains the info for descending at min_class
    end;
    erase_uniques( uniques, nonzerocol, candidate, images );
    if impossible then return possibilities; fi;
    
    # step 3: collapse
    
    remain:= Filtered( [ 1 .. nccl ], x -> nonzerocol[x] );
    for i in [ 1 .. Length( matrix ) ] do
      matrix[i]:= Sublist( matrix[i], remain );
    od;
    candidate:=     Sublist( candidate, remain );
    nonzerocol:=    Sublist( nonzerocol, remain );
    moduls:=        Sublist( moduls, remain );
    matrix:= ModGauss( matrix, moduls );
    ncha:= Length( matrix );
    pos:= 1;
    fusionperm:= [];
    newimages:= [];
    for i in remain do
      fusionperm[ i ]:= pos;
      if IsBound( images[i] ) then newimages[ pos ]:= images[i]; fi;
      pos:= pos + 1;
    od;
    min_class:= fusionperm[ min_class ];
    for i in Difference( [ 1 .. nccl ], remain ) do
      fusionperm[i]:= pos;
      newimages[ pos ]:= images[i];
      pos:= pos + 1;
    od;  
    images:= newimages;
    fusion:= CompositionMaps( fusionperm, fusion );
    nccl:= Length( nonzerocol );
    
    # step 4: backtrack
    
    evaluate:= function( candidate, nonzerocol, uniques, images )
    local i, j, col, val, row, quot, abgespalten, step, erster,
    descendclass, oldimages;
    abgespalten:= erase_uniques( [ uniques ], nonzerocol, candidate, images );
    if impossible then return abgespalten; fi;
    descendclass:= min_class;
    oldimages:= images[ descendclass ];
    for i in [ 1 .. min_anzahl ] do
      images[ descendclass ]:= oldimages[i];
      oldrows:= evaluate( Copy(candidate), Copy(nonzerocol), descendclass,
                          Copy(images) );
      Append( matrix, oldrows );
      if Length( matrix ) > ( 3 * ncha ) / 2 then
        newmatrix:= [];            # matrix:= ModGauss( matrix, moduls );
        for j in [ 1 .. Length( matrix[1] ) ] do
          if nonzerocol[j] then
            row:= StepModGauss( matrix, moduls, nonzerocol, j );
            if row <> false then Add( newmatrix, row ); fi;
          fi;
        od;
        matrix:= newmatrix;
      fi;
    od;
    return abgespalten;
    end;

    descendclass:= min_class;
    oldimages:= images[ descendclass ];
    for i in [ 1 .. min_anzahl ] do
      images[ descendclass ]:= oldimages[i];
      oldrows:= evaluate( Copy(candidate), Copy(nonzerocol), descendclass,
                          Copy(images) );
      Append( matrix, oldrows );
      if Length( matrix ) > ( 3 * ncha ) / 2 then
        newmatrix:= [];             # matrix:= ModGauss( matrix, moduls );
        for j in [ 1 .. Length( matrix[1] ) ] do
          if nonzerocol[j] then
            row:= StepModGauss( matrix, moduls, nonzerocol, j );
            if row <> false then Add( newmatrix, row ); fi;
          fi;
        od;
        matrix:= newmatrix;
      fi;
    od;
    return possibilities;
    end;


#############################################################################
##
#F  ContainedCharacters( <tbl>, <constituents>, <paracharacter> )
##
##  the list of all characters in <paracharacter> which lie in the linear span
##  of the rational characters <constituents> (modulo centralizer orders) and
##  have nonegative scalar products with all elements of <constituents>.
##
##  (The elements of the returned list are not necessary linear combinations
##  of <constituents>.)
##
ContainedCharacters := function( tbl, constituents, paracharacter )
    local func, chi, choice;

    func:= function( chi )
    return NonnegIntScalarProducts( tbl, constituents, chi ); end;

    if IsInt( paracharacter[1] ) then
      choice:= [];
      for chi in constituents do
        if chi[1] <= paracharacter[1] then AddSet( choice, chi ); fi;
      od;
      constituents:= choice;
    fi;
    return ContainedDecomposables( constituents, tbl.centralizers,
                                   paracharacter, func );
    end;


